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Abstract. The impact of various competing heating pro- 
cesses on the thermal evolution of neutron stars is investi- 
gated. We show that internal heating leads to significantly 
enhanced surface temperatures for pulsars of middle and 
old age. The heating due to thermal creep of pinned vor- 
tices and due to outward motion of proton vortices in the 
interior of the star leads to a better agreement with the 
observed data in the case of enhanced cooling. The strong 
pinning models are ruled out by a comparison with the 
cooling data on the old pulsars. For millisecond pulsars, 
the heating due to thermal creep of pinned vortices and 
chemical heating of the core have the largest impact on 
the surface temperatures. The angular dependence of the 
heating rates require two dimensional cooling simulations 
in general. Such a simulation is performed for a selected 
case in order to check the applicability of one-dimensional 
codes used in the past. 
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1. Introduction 

Simulations of the thermal evolution of neutron stars con- 
fronted with the soft X-ray and extreme UV observations 
of the photon flux emitted from their surface provide one 
of the most useful diagnostics of the physical processes 
operating in the dense interior of such stars. The early 
evolution of a neutron star is completely dominated by 
the cooling via neutrino emission; at this stage the star 
effectively radiates away the initial content of its thermal 
energy. In the later epoch one faces a competition between 
cooling and heating processes, where the heating may even 
dominate in the late stages of the thermal evolution. Even- 
tually, they reach a thermal equilibrium phase, where the 
heat generated in the star is radiated away at the same 
rate from the stellar surface. 



The relative importance of the various elementary 
processes contributing to cooling is fairly well un- 
derstood (e.g. Richardson et al., 1982; Van Riper, 1991; 
Umeda et al., 1994; Schaab et al., 1996; Page, 1998). The 
investigations performed so far have however the deficit 
that they did not sufficiently clarify the relative weight of 
the different heating processes. The majority of these in- 
vestigations concentrated on few preferred dissipation pro- 
cesses (e.g Shibazaki & Lamb, 1989; Cheng et al., 1992; 
Sedrakian & Sedrakian, 1993; Reisenegger, 1995) with 
the exception of the work of Van Riper (1991), who treats 
the heating in general phenomenological terms. A mean- 
ingful comparison of the individual processes studied in 
these papers is strongly hampered by the fact that these 
investigations were performed for different microphysical 
inputs and different levels of sophistication concerning 
the cooling simulation code, preventing one from draw- 
ing definitive conclusions about the relative importance 
of the individual heating processes. 

The aim of this paper is a comparative analysis of 
the impact of different heating processes on the thermal 
evolution of neutron stars. Our simulations take into ac- 
count the most recent developments concerning the mi- 
crophysical input and employ both relativistic and non- 
relativistic equations of state of superdense matter. We 
perform the cooling calculations for a fully relativistic, 
evolutionary simulation code. Although internal heating 
is generally azimuthally asymmetric, we shall use an one- 
dimensional cooling code for almost all simulations. To 
demonstrate the reliability of this approximation, a two- 
dimensional calculation is performed for a selected case as 
well (Schaab & Weigel, 1998). As far as we know, this is 
the first fully two-dimensional simulation of neutron star 
cooling which takes internal heating into account. 

The heating processes emerge as a response to 
the loss of rotational energy of the star and can be 
divided roughly into two groups. Firstly, processes 
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cnist/corc-quakcs (Cheng et al., 1992 and references 
therein); nuclear reactions in cold nuclear matter with 
non-equilibrium composition (picno-nuclear reactions 
in the crust and weak processes in the core, see, e.g., 
lida & Sato, 1997 and Reisenegger & Goldreich, 1992). 
Secondly, processes related to superfluidity such as 
dissipativc motion of the neutron vortex lattice 
in the superfluid crusts (Shibazaki & Lamb, 1989; 
Link & Epstein, 1991; Bildsten & Epstein, 1989; 

Jones, 1990) and the dissipativc motion of proton 
vortex lattice in the superfluid core and the vortic- 
ity decay at the phase separation boundaries, e.g. at 
the crust-corc interface, (Sedrakian & Sedrakian, 1993; 
Sedrakian & Cordes, 1998b). 

The rate of dissipation associated with these processes 
is, in most cases, intimately connected with the spin evolu- 
tion of the star, in particular the evolution of its magnetic 
moment. The time dependence of the magnetic moment 
and the moment of inertia will affect the heating rate and 
thus the star's cooling behaviour along with the relation 
between the luminosity and the age of the star; these ef- 
fects are ignored in this study in order to constrain the 
actual number of the evolutionary scenarios; we shall as- 
sume therefore a constant magnetic moment and moment 
of inertia. 

The paper is organised as follows. The equations of 
thermal and rotational evolution of a neutron star are 
given in Sect. 2. Various heating processes associated with 
the dissipative motion of vortices and the changes of the 
equilibrium structure are reviewed in Sects. 3 and 4. Addi- 
tional physical input is discussed in Sect. 5. The observed 
data are summarised in Sect. 6. Section 7 presents the 
results of our cooling simulations which are discussed in 
Sect. 8. 



2. Thermal and rotational evolutions 

2.1. Equations of thermal evolution 

The thermal evolution is governed by the coupled system 
of equations for energy balance (Thorne, 1977), 



a(Le2'^) 



dr 



and thermal energy transport. 



dr 



4nr'^K 



(2) 



This system requires as a microphysical input the neutrino 

cmissivity per unit volmnc, €,j(p, T), the heat rate per unit 
volume, h(p,T,n,^), the heat capacity per unit volume, 
Cv(p,T), and the thermal conductivity, K,(p,T), where p 



e , respectively, L is the luminosity, and the other vari- 
ables have their usual meaning. The boundary conditions 
for (1) and (2) read 



L(r = 0) = , 



Tm{rm, Lra, M^) , 



(3) 
(4) 



where T^irm, L^, M^) is fixed by the properties of the 
photosphere at the density p = pm = 10^°gcm~^ 
(Gudmundson et al., 1983). Since the star's thermal evo- 
lution for times greater than, say, a few months does not 
depend on the detailed initial temperature profile in the 
interior of the star, one can choose, without loss of gen- 
erality, an initialy isothermal temperature distribution of 
T(r) = 10^^ K. The set of equations for thermal evolution 
(1) and (2) were solved numerically by means of a Ncwton- 
Raphson type algorithm. The microphysical input quan- 
tities needed for the solution of these equations are given 
in Table 1 and will be discussed in detail in Sect. 5. If one 
accounts for azimuthal asymmetry of the input quanti- 
ties, e.g. the heat rate h, one has to solve two-dimensional 
thermal evolution equations, which have already been pre- 
sented in Schaab & Weigel (1998) Schaab98c along with 
the numerical method for solving them. 

2.2. The normal component 

We shall start the discussion of the rotational evolution by 
a brief review of the spin dynamics of the normal (non- 
superfluid) component of the star. The combined effect of 
the braking of the star's rotation via the emission of elec- 
tromagnetic radiation, electron-positron wind, and gravi- 
tational waves on the spin evolution is commonly cast in 
a generic power-law form 



Oc(t) = -K{t)n^{t), 



(5) 



where flc is the spin frequency of the normal compo- 
nent. The fimctional form of K{t) and n{t) depends on 
the dominant process assumed. If one assumes that the 
energy losses are due to the magnetic dipole radiation 
(Goldreich et al., 1971), then n = 3 and is a fimction of 
the surface magnetic field strength B, the angle a between 
spin and magnetic axis, the moment of inertia of the star 
I, and a correction factor 7 = 1 -|- dln7/dlnf2c: 



K(x 



B"^ siv? a 



(6) 



The braking index is known for the four youngest pulsars 
for which the measurements of the second time deriva- 
tive of the period, P, are available (Lyne et al., 1993, 
Kaspi et al., 1994, Boyd et al., 1995). All four braking 
indices given by n = PP/P"^ fall into the range 
between 1.2 and 2.8. The deviations from the value 
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(Blandford et al., 1983, Muslimov & Page, 1996), an in- 
crease of the angle between the spin and magnetic mo- 
ment axis (Beskin et al., 1984, Link & Epstein, 1997), or 
presence of the weakly coupled superfluid in the core of 
the star (Sedrakian & Cordes, 1998a). We shall assume, 
in a first approach to the problem, K to be constant and 
n = 3. 

With these assumptions, integration of Eq. (5) yields 

(n-l)i^i + fir("-')) , (7) 

where Cli is the initial angular velocity at f = 0. If the 

initial angular velocity is large compared to the present 
angular velocity, Eq. (7) simplifies to 

OeW = ((n- (8) 
and the age of the pulsar can be expressed as 

t= . ^ . (9) 

n-i\n,{t)\ 

Apart from affecting the cooling history, the time varia- 
tion of K, as well as a deviation of the braking index from 
its canonical value, n = 3, may lead to over- or underesti- 
mating the true age by a factor as large as 3. A compari- 
son of the theoretical cooling tracks with observed data is 
therefore affected by this uncertainty (see Sect. 6). 

2.3. The superfluid component 

The angular velocity of the superfluid phases adjusts to 
the changes in the rotation rate of the normal compo- 
nent via expansion or contraction of the neutron vortex 
lattice, which carries the angular momentum of the su- 
perfluid. The mean macroscopic density of neutron vortex 
lines ni in cylindrical coordinates (rp, z, ({)) is given by the 
Feynman-Onsager quantisation condition 

Km{rp,t) = 2a(rp,i) +rp^^^^^. (10) 

f7s(? p) denotes the angular velocity of the superfluid com- 
ponent and K — h/2m^ is the quantum of circulation 
carried by each vortex, where h denotes the Planck con- 
stant and rriN is the bare neutron mass. The second term 
accounts for the gradients in the distribution of vortex 
lines. Because of the lack of sources and sinks for vortex 
lines in the bulk of the superfluid, the temporal evolution 
of the vortex number density obeys the number conserva- 
tion law: 

^ + V-{nivi) = Q, (11) 

where vi is the velocity of vortex lines. Combination of 
Eqs. (10) and (11) leads to 



A decrease of the spin of the superfluid component {(ig < 
0) causes an outward radial velocity i)[ > 0. In the follow- 
ing we neglect the second term in parentheses and assume 
r2s(rp) to be constant throughout the superfluid phase, 
since the vorticity gradients depend on the density pro- 
file of the mutual friction coefficients, which are poorly 
known. Furtlicairiorc, over the evolutionary timcscalc one 
can assume that tic = tls, i.e. the deceleration rates of the 
superfluid and the normal components are the same. Then 
the velocity of the radial motion of the vortices (and there- 
fore the rate of dissipation) are related to the spin-down 
characteristics of the normal component via Eq. (12). In 
this manner the dynamics of superfluid phases is elimi- 
nated from the system of coupled equations of secular and 
thermal evolution, and the heating rates are given directly 
in terms the spin characteristics of the normal component. 

3. Dissipative motion of vortex lattices 

We next turn to the description of the dissipation pro- 
cesses included in the thermal evolution simulations. A 
number of dynamical coupling regimes has been suggested 
in the literature for the coupling of the superfluid compo- 
nent to the crusts. We shall discuss several variants of 
the two main alternatives - the pinned and the corotation 
regimes. 

3.1. Vortex pinning at nuclei 

The vortex creep theories of the dynamics of neutron 
vortices in the inner crust assume these vortices to be 
pinned to the nuclear lattice (Anderson & Itoh, 1975, 
Alpar, 1977, Alpar et al., 1984, Epstein & Baym, 1988, 
Link & Epstein, 1991, Pizzochero et al., 1997). The pin- 
ning force provides a large resistive force to the vortex 
motion which is necessary to spin down the superfluid 
component. In this regime, the vortices move due to the 
process of vortex creep, i.e. by thermal activation from 
one pinning configuration to another. The vortex creep 
rate has the general form (Link et al., 1993) 

vl = voe-^/^'^' , (13) 

where vq is a microscopic velocity, A denotes the activation 
energy for a segment of a vortex line to overcome its pin- 
ning barrier, and Teg is the effective temperature which 
includes quantum tunnelling effects and is a function of 
the local temperature of the thermal bath. These quanti- 
ties depend on microscopic parameters (e.g. the pinning 
energy), which themselves depend on the density, and on 
the velocity difference between the superfluid and the nor- 
mal component 5v'^. Since the creep rate v\ is given by Eq. 
(12), Eq. (13) has to be inverted to obtain the velocity dif- 
ference 5v'^. This is done by an iterative procedure similar 
to the one described by Van Riper et al. (1995) . 
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The energy dissipated by a single vortex per unit time is 
-Ediss = Jm'^i, where = PskSv"^ is the radial component 
of the Magnus force. The dissipated energy h per unit time 
and volume is then given by: 



h = E'dissm = TpPsSv'^lCta 



(14) 



where Eq. (12) has been used. The integration over the 
volume of the superfluid component gives: 



J hdv = J , 



arp5v'''dV. 



(15) 



This result guarantees eonscrvation of both energy and 
angular momentum. For the implementation of h in our 
one-dimensional cooling code, h has to be averaged over 
spherical shells: 



h 



Jhsmed0d(j) _ 1 



Air 



-nrpsSv'''\Qs\, 



(16) 



where we assumed that Sv"^ is approximately constant 
on spherical shells. This approximation seems to be re- 
liable because Sv'^ depends only logarithmically on rp 
(Van Riper ct al., 1995). 

We shall consider two models for pinning: the first 
model (EB) uses the parameters from Epstein & Baym 
(1988) and Link & Epstein (1991) , where the pinning 
energy is derived from the Ginzburg-Landau theory. In 
the second model (PVB) the parameters are derived 
by Pizzochero et al. (1997) from the Bardeen-Cooper- 
Schrieffer theory of superconductivity. Below the den- 
sity 1.5 X lO^^gcm"^ the parameters for the both mod- 
els do not deviate considerably. However, above densities 
> 1.5 X lO^^gcm""^ the large deviations in the pinning 
energy (the typical pinning energy is ~ 10 MeV in the EB 
model and ^ 1 MeV in the PVB model) leads to rather 
different evolutionary scenarios and, therefore, we shall 
treat both models seperately. 



3.2. Corotating regime in the crust 

The dynamics of the neutron vortices in the corotating 
regime is governed by the force balance equation 



/m + fd 



0, 



(17) 



where the Magnus and the drag (resistive) forces are: 

fu = -psK' X {Vs - Vl), (18) 

fd = -v{vi-v,), (19) 

with the drag coeflacient r]. A decomposition of the force 
balance equation into radial and azimuthal components 
gives 



where the dissipation angle ^d) following Bildsten & Ep- 
stein (1989), is defined as 



tan^d = . 

PsK 



(21) 



The rate of the dissipation per unit volume is then given 

by 



h = - fd ■ {Vl - Vc)nfree = Ps^-plfis 



CO, 



with the abbreviation 



Vs - Vc 



Tp sin 6d cos 9d ' 
1 1 



2f2s $free sin 6d cos ^d 



(22) 

(23) 
(24) 



Here, $free denotes the fraction of the free (corotating) 
vortices. This quantity is largely unknown, but for sta- 
tionary situations, one may assume that definite crustal 
regions support either corotation or creep regimes and the 
mixing is small. Hence we shall adopt $free '■= ?^free/'^i = 1 
for the corotating regime. 

The dominant processes in this regime result from the 
coupling of the translational motion of the vortices to the 
electron-phonon system (Jones, 1990). The leading order 
coupling implies a resistive force constant 



aE?. 



(25) 



where a is the lattice constant, Ep the pinning energy, 
£,n the neutron coherence length, M the (effective) mass 
of the nuclei, and Cg the phonon velocity (given by the 
unperturbed dispersion relation). Using eq. (25) we find 
rjo ~ 10^ g cm^^ s^^ at the density Ps = 10^^ g cm~^ 
and, therefore, tanOd ~ 5 x 10~^. The processes cor- 
responding to higher orders in the vortex displacement 
are generally velocity dependent. The next-to- leading or- 
der term, involving one-phonon processes, is estimated as 
(Epstein & Baym, 1992, Jones, 1992) 



El 



2psu;a^l V27r|nC^ 



(26) 



where ck is the kelvon velocity due to excitation of the 
oscillatory degrees of freedom of the vortex lines moving 
with velocity v^. This process is sensitive to the actual ex- 
citation spectrum of kelvons. The possible localisation of 
kelvons in a random nuclear potential (Jones, 1992) leads 
to a threshold in the excitation spectrum beyond which 
a continuum of states becomes available for excitations. 
In addition one requires knowledge of the scattering am- 
plitudes beyond the Born approximation. For typical pa- 
rameters one finds rji ~ 10"'^'^ g cm"-'^ s~^ at a density of 
Ps ^ 10^'^ g cm~'^, and, therefore, tan^^ ~ 1. However, in 
the stationary situations, the vortex velocities are rather 



Schaab et al.: Impact of internal heating on the thermal evolution 



5 



3.3. Electron-vortex scattering in superfluid core 

An estimate of the heating processes in the superfluid core 
requires knowledge of vortex lattice structure that nucle- 
ates after the superfluid phase transition. We shall pur- 
sue the point of view that the core's magnetic field nu- 
cleates in a proton vortex lattice. The frictional forces 
in a neutron star's core are dominated by the scatter- 
ing of the quasiparticle excitations of the electron nor- 
mal Fermi-liquid by the proton vortex lattice. The rea- 
son is that the number of proton vortices per neutron 
vortex, given by n ~ (\/3/2)(B/$o)c?ra, is of the order 
of n > 10^^ for B ~ 10^^ G and a neutron vortex lat- 
tice constant of dn ~ 10~^ cm. The dominant dissipa- 
tive scattering results from the scattering of electron ex- 
citations off the quasiparticle excitations confined in the 
core of proton vortex lines. The scattering rate is a sen- 
sitive fimction of the form-factor of the vortex line core 
and the actual distribution of vortex lines (e.g. formation 
of vortex clusters, Sedrakian & Cordes, 1998b). An order- 
of-magnitude estimate of the dissipation rate is given by 
(Sedrakian & Sedrakian, 1993) 

h^^fT(^] R^sm^6cos^e, (27) 

where T is the temperature, R the radius of the superfluid 
core and / is defined as^. 



/ = 4l7rV"'P 



h^C fcpT 



2mpC^ 



1/2 



exp 



2-|-2/3|fe| 



0.78A2 

EpkBT 



(28) 



where fcp — (37r^rip)^/'* denotes the proton Fermi momen- 
tum, kpT = {4:kFmpe^ /-Kh^Y^"^ the Thomas- Fermi screen- 
ing length, Ep = fi,^fcp/2m* the proton quasiparticle en- 
ergy with the proton effective mass m*, k = (m*— mp)/mp 
the entrainment coefficient, and 



h^kF 
TrmpAp ' 



Ap = (l + |fc|)^/' 



JTlpC 



1/2 



(29) 



(30) 



the proton coherence length and the magnetic field pen- 
etration length, respectively. It should be noted that the 
dynamical coupling of the proton vortex clusters to the 
electron liquid is controlled by the electron flux-scattering 
processes, whose rate is larger than that for the electron 
scattering off the vortex core quasiparticles. While in the 
latter case the scattering is inelastic (i.e. energy is trans- 
fered to the proton quasiparticle localised in the vortex 

^ Note that the multiDlication sign in the second line of Ea. 



core), the former process of electron- flux scattering is pre- 
dominantly elastic i.e. involves only a change of direction 
of the electron momentum (but no changes in the energy). 
In the classical (or quasi-classical) limit this process cor- 
responds to the bending of the electron trajectory in the 
magnetic fleld of a vortex under the action of the Lorentz 
force and is manifestly non-dissipative. In the quantum 
limit the processes of the soft-photon emission can not 
be dismissed in general; however the corresponding rates 
should be substantially less than the elastric scattering 
rates. 

3.4. Decay of vortices 

The proton vortices, which are dragged along by the neu- 
tron vortex motion in the radial outward direction to the 
crust-core boundary, would decay as they merge in the 
crust, thus releasing their self-energy. The relevant crust- 
core boundary for this process should be identified with 
the phase transition point between the phase with pro- 
tons bound in heavy nuclei and the phase in which they 
arc in contimmm states. Since the crustal matter can not 
support a proton supercurrent, the circulation currents 
of proton vortices would decay on the ohmic dissipation 
timcscalc, which is much smaller than a star's spin-down 
time scale. The resulting heating rate at the crust-core 
boundary is given by 



/ 



hdV 



with 




(31) 



2/3|fe| 



V47rAp 



In^ 

SP 



(32) 



Here, $0 = hc/2e is the flux quantum. The quantity g is 
very sensitive to variations of the microphysical parame- 
ters such as the proton density, the effective mass, and the 
gap energy at the location of the phase transition point. 
We use a fixed upper limit for this process corresponding 
to the value of 3 = 10^^. 

4. Readjustment to equilibrium structure 

As a neutron star spins down its global structure readjusts 
itself toward a more spherically symmetric configuration 
to minimise the sum of gravitational and rotational en- 
ergy. Two models of internal heating discussed here are 
based on the structural readjustment of the oblateness of 
the star, defined as 
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Fig. 1. Oblateness e and the variation of the central density 1 — 
nc(n)/nc(0) of rotating neutron stars versus angular velocity. 
The parameters are: Mb = 1.58Mq and UVm+UVII-EOS. 

and of the central baryon density ric- /(O) denotes the 
star's moment of inertia. The equihbrium values of both 
quantities arc determined as functions of the angular 
velocity by computing a sequence of rotating neutron 
star models with constant baryon number (Schaab, 1999). 
These calculations arc done by solving the general rcla- 
tivistic stellar structure equations in two dimensions, in a 
similar manner as in Komatsu et al. (1989) and Bonaz- 
zola et al. (1993). As an example the oblateness e and the 
fraction 1 — nc(n)/nc(0) with the central density nc{^) 
are plotted for a neutron star model with constant baryon 
mass Mb = 1.58Mq in Fig. 1. The underlying equation of 
state is the non-relativistic UV14+UVII- model (see Sect. 
5.1). 

4-1- Crust cracking 

If both the interior and the crust are liquid, the oblate- 
ness will continuously adjusts itself to the equilibrium 
value. This is the case for the temperatures T > 10^ K 
(Cheng et al., 1992). However, if the crust solidifies, the 
oblateness continuously depart from its equilibrium value 
causing an increasing mean stress on the crust according 
to (Baym & Pines, 1971) 



a = fl{e- Ceq), 



(34) 



where jl is the mean shear modulus of the crust. The quan- 
tity Ceq refers to the equilibrium value of the oblateness. 
If the mean stress a reaches some critical value 



(35) 



the crust breaks and the deviation of e from its equilibrium 
value is reduced. For an ideal Coulomb lattice, the value 



of impurities or defects in the neutron star crusts could 
be considerable due to, e.g., the fast cooling below the 
melting temperature of the lattice. Smolukowski (1970) 
argues that the value of 6c is then reduced to roughly 
10"'' - 10"^ 

During the breaking of the crust the strain energy 



2B{e 



)A(e 



(36) 



is released (Cheng et al., 1992), where A(e — eeq) <C (e — 
Ceq) is the change in (e— Ceq), and the coefficient B is given 
in terms of the mean shear modulus by B = ■^Vcfi, with 
Vc being the volume of the crust. We follow Cheng et al. 
(1992) in assuming that (e — Ceq) is approximately equal 
to its maximum value 9c and that the time At between 
two successive quakes is small compared to the time scale 
of thermal evolution. Hence^ 



dK, 



strain 



dt 



At 



2B0, 



A(e - eeq) 

At 



dCcq 
-zBuc~r^ — i^r 



d^c 



(37) 



The crust strength B is given by the integral 
(Baym & Pines, 1971) 



-i 



B= / 6(r)dV 



(38) 



with 



^(r) = ^,C,4r)(96-U6(L)\IIlQ^ (39) 



and the shear modulus of a bcc-lattice 
(Mott & Jones, 1936) 



Cii{r) = 0.3711Z^e^nN 



-1/3 



(40) 



The total strain energy release (37) can thus be expressed 
in terms of a local energy loss rate: 



h = 2b0c^\nc\. 



dflc 



(41) 



4-2. Chemical heating 



The central density of a rotating neutron star with fixed 
baryon number depends on the angular velocity (see Fig. 
1). As the star spins down, the centrifugal force decreases 
and correspondingly the central density increases (by up 
to 20 %). This is accompanied by a shift of the chemical 
equilibrium. The matter would maintain chemical equilib- 
rium if the relaxation timescales for the weak processes are 
small compared to the timescale of rotational evolution. 

^ Our formulas differ from those of Cheng et al. (1992) by a 

factor 6c/ e, which is <C 1 for slowly rotating stars. We there- 
fore obtain no effect on the coolinE of slowlv rotatine stars, in 
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However, these timescales are found to be comparable 
and the actual composition departs from chemical equi- 
librium, which modifies the reaction rates (Haensel, 1992) 
and usually leads to a net conversion of chemical en- 
ergy into thermal energy (Reisenegger & Goldreich, 1992, 
lida & Sato, 1997). 

For simplicity wc restrict ourselves to a system con- 
taining only neutrons, protons, electrons, and muons in 
/3-equilibrium (as is described by the non-relativistic EOS 
UVi4H-UVII). The equation of state of cold, charge neutral 
matter, which is not necessarily in /3-equilibrium (we as- 
sume however equal chemical potentials for electrons and 
muons: /ic = l^tJ.-), therefore depends only on two param- 
eters. These are chosen to be the baryon density n and 
the relative proton density x = rip/n. These two param- 
eters are constant on closed, sphere-like surfaces S{N) of 
constant energy density p and pressure P. N refers to the 
number of baryons enclosed by S. We approximate the 
time derivative of the baryon density on a specific surface 
S{N) by the expression (Reisenegger & Goldreich, 1992) 



dn 
dt 



n 

Tie 



dric di7 
djT d^' 



(42) 



where ric is the central baryon density, and dric/ dO is ob- 
tained from calculations of rotating neutron star sequences 
(see Fig. 1). 

The equilibrium value a;eq of the proton fraction is 
determined by the following condition on the energy per 

baryon E: 



5fi := 



dE dE dup dE dn^ dE dn^ dE dn 



dx 



drip dx 



+ 



drin dx dne dx 



+ 



dn^ dx 



Mp - Mn + Me = 0, 



(43) 
(44) 



where we used the charge neutrality of the system, 
dup/dx = dn^/dx + dn^/dx, the condition jif. = ji^, and 
the equation dup/dx = —dnn/dx = n. The time deriva- 
tive of n is linked to the time derivative of x'^'^ by 



dn 



0=0 



dn dO 
di^ dF' 



(45) 



where Eq. (42) was used. 

A deviation from chemical equilibrium corresponds to 

a non- vanishing value of 5^. The non- vanishing gradient 
of the total chemical potential including both internal and 
external contributions causes a diffusion of particles. If 
the timescale for chemical relaxation is larger than the 
timescale for this diffusion, the chemical relaxation takes 
place in a larger regoin of the neutron star, and (5/i has 
to be averaged over such regions. It seems likely that 
this kind of relaxation occurs in the core of neutron stars 
(Reisenegger, 1997). 

The time derivative of 



yields the evolution equation for Sji: 
d'^E 



dSp, 
'df 



dx^ 



da;'*'' n 
dn rir 



dn df2 
dO 



(47) 



where r(^;Lt, n) = ndx/dt is the difference between the 
rates per unit volume of electron capture and beta de- 
cays. This difference gives a positive value for 5^ > 0, so 
that chemical equilibrium tends to be restored. The sec- 
ond derivative of the energy per baryon is equal to 



d^E 



dx"^ 



= 85(n) + \hc{^TmY'^x-^'^, (48) 



where S{n) denotes the symmetry energy. 

The average energy released by these reactions into 
kinetic and thus thermal energy is equal to the difference 
of the Fermi energies of protons and electrons on the one 
side and neutrons on the other side, i.e.: 



h = T5n. 



(49) 



At the same time, thermal energy is lost by emission of 
thermal neutrinos, as it is the case for matter in chem- 
ical equilibrium. The emission rates are however mod- 
ified by the deviation from chemical equilibrium (see 
Reisenegger & Goldreich, 1992 for the values of F and Ci,)- 
Like the neutrino emissivity Cj^, the reaction rate F is sup- 
pressed by superfluidity unless the deviation of chemical 
equilibrium Sp, exceeds the sum of the gap energies of the 
participating baryons (Reisenegger. 1997). 

Until now we have considered only the restoration of 
chemical equilibrium in the core. The situation in the crust 
is similar, though more complex. lida & Sato (1997) stud- 
ied the Lagrangian changes in pressure associated with ele- 
ments of matter due to the spin-down in the framework of 
the Hartle-approximation (Hartle, 1967). By considering 
the nuclear processes induced above pi = 3 x 10^^ g cm~'^ 
by compression or decompression of matter, they approx- 
imate the chemical energy per unit volume and time con- 
verted to thermal energy to 



hKnqNsm.e\l-lA&cos^e\ 



2n\VL\ 



(6283s-i)2 



(50) 



where nq « 4eVfm~^ is the energy per unit volume 

released by one non-equilibrium process, N « 15(p — 
Pi)/(ptr — Pi) is the number of processes in a unit cell, 
and ptr = 1-7 x lO^^gcm"^ is the transition density be- 
tween core and inner crust. By averaging over spherical 
shells we obtain 
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5. Additional ingredients and competrison of the 
efficiency of the heating mechanisms 

Apart from the rates of the various heating processes dis- 
cussed in the previous sections, some further ingredients 
are needed to solve the general relativistic equations of 
stellar structure and evolution (see Sect. 2.1). These addi- 
tional ingredients are summarised in Table 1. Here we shall 
only discuss the equation of state of the dense interior of 
a neutron star, the neutrino emissivities, the superfluid- 
ity gaps and the photosphere, since these ingredients are 
the most important ones. Further, we will compare the 
efficiencies of the various heating mechanism. 

5.1. Equation of state 

For the outer and inner crust wc adopt the equations 
of state of Haensel & Pichon (1994) and Negele & Vau- 
therin (1973). The transition density between the ionic 
crust and the core of a neutron star is taken to be ptr = 
1.7 X 10^'*gcm-3 (Pethick et al., 1995). For the present 
study, we choose two models of high density matter. The 
first, non-relativistic model (UV14+UVII) is obtained by 
solving the Schrodinger equation by means of a variational 
approach (Whinga et al., 1988). The other model (RHF8) 
uses relativistic Briickner-Hartree-Fock results up to 2- 
3 times normal nuclear density. Hyperons are included 
within the relativistic Hartree-Fock approach above this 
density (Huber et al., 1998). One important difference be- 
tween these equations of state is that the non-relativistic 
model treat neutron star matter as being composed of 
neutrons and protons only (which are in /3-equilibrium 
with leptons), whereas the relativistic model accounts for 
all hyperon states that become populated in the cores of 
neutron stars. 

5.2. Neutrino emissivity 

The neutrino emission processes can be divided into slow 
and enhanced processes depending on whether one or two 
baryons participate in the reaction. Due to the rather dif- 
ferent phase spaces associated with both kind of processes 
the emission rates differ by several orders of magnitude. 

If one neglects exotic states in neutron star matter (like 
quark-gluon plasma and meson condensates) the only en- 
hanced neutrino emission processes are the direct nucleon 
and hyperon Urea processes 

n^p + \- +i?x (52) 
Bi^B2 + X-+ux, (53) 

respectively, where _Bi.2 = n, p. S^'*^, A, S"'^ denotes the 
baryons and Ai,2 = e~, /x~ the leptons. Because of the 
/3-equilibrium the inverse reaction (with D replaced by i^) 



Simultaneous conservation of energy and momentum re- 
quires that the triangle inequality < p^^ and the 
two inequalities obtained by cyclic permutation are ful- 
filled for the Fermi momenta pj . If the inequalities are not 
fulfilled the process is extremely unlikely to occur and the 
corresponding emissivity can be neglected. The availabil- 
ity of the various fast processes thus depends on the partial 
concentrations of each baryon species. The proton frac- 
tion, which determines the possibility of the direct nucleon 
Urea process, depends crucially on the symmetry energy, 
which is poorly known for high density matter. As an ex- 
ample of an equation of state which yields slow cooling we 
study tlie non-relativistic equation of state UV14+UVII. 
Due to the non-monotonic behaviour of the symmetry en- 
ergy at high densities, the triangle inequality is not ful- 
filled for this eriuation of state. The relativistic equation of 
state RHF8 allows for both the direct nucleon and the di- 
rect hyperon Urea processes for star masses M > 1.09 Mq 
and M > 1.22 Mq, respectively. 

During the transition to the superfluid state of neu- 
trons and protons the superfluid pair breaking and forma- 
tion processes become important (Flowers & Itoh, 1976, 
Voskresenskii & Senatorov, 1987). We refer to Schaab et 
al. (1997) for a detailed description of these processes and 
their impact on the cooling of neutron stars. 

5.3. Superfluidity 

The pairing gaps are sensitive to the underlying micro- 
scopic model for the nucleon-nucleon interaction. Since 
these models are rather uncertain, especially at high den- 
sities, the value of the maximum gap and the density range 
where pairing can occur is model dependent. The disagree- 
ment arising from the different interaction models is en- 
larged by more subtle issues of the medium polarisation 
(Wambach et al., 1991, Schulze et al., 1996) or the inclu- 
sion of relativity (Elgar0y et al., 1996). 

At densities of the order 2 x lO^^gcm"^ the S 
state interaction becomes repulsive for neutrons, and 
the ^So neutron gap closes. The attractive — ^ F2 
state interaction leads to the pairing in this channel at 
about the same densities. For the P state pairing we 
use the recent calculation of Baldo et al. (1998) which 
is based on the modern nucleon-nucleon CD-Bonn po- 
tential (Machleidt et al., 1996), which provides accurate 
fits to the actual nucleon-nucleon scattering data below 
350 MeV. 

Field theoretical descriptions of the heavy-baryon in- 
teraction based on the one-boson-exchange model indicate 
that the A — A interaction has an attractive component. 
From the so-called doubly-strange hypernuclei, where two 
hyperons are boimd in a single nucleus (Imai, 1992), one 
can verify this field theoretical result, since, for instance, 
the separation energy of the two A hyperons. Baa, ex- 
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Parameter 


References 


Equations of state: 




crust 


Haensel & Pichon, 1994, Negele & Vautherin, 1973 


core (alternatives) 




UV14+UVII 


Wiringa et al., 1988 


RHF8 


Huber et al., 1998 


Superfluidity 


see Table 2 


Heat capacity 


Shapiro & Teukolsky, 1983, Van Riper, 1991 


Thermal conductivities: 




crust 


Itoh et al., 1984, Itoh, 1983, Mitake et al, 1984 


core 


Gnedin & Yakovlev, 1995 


Neutrino emissivities: 




pair-, photon-, plasma-processes 


Itoh et al., 1989 


bremsstrahlung in the crust 


Yakovlev & Kaminker, 1996, Haensel et al., 1996 


bremsstrahlung in the core 


Priman & Maxwell, 1979, Kaminker et al., 1997 


modified Urea 


Friman & Maxwell, 1979, 




Yakovlev & Levenfisli, 1995 


direct nucleon Urea 


Lattimer et al., 1991 


direct hyperon Urea 


Prakash et al., 1992 


superfluid pair breaking 


Voskresenskii & Senatorov, 1987, 




Schaab et al., 1997 


Photosphere: 


Potekhin et al., 1997 



Table 2. Gap energies of superfluid states in neutron star matter. The density ranges are calculated for the RHF8 equation of 
state. 



Pairing 


Amax [MeV] 


Density range [fm 


References 


neutron ^So 


1.01 


< 0.16 


Sehulzc ct al., 1996 


proton ^So 


0.92 


0.10-0.31 


Elgar0y et al., 1996 


neutron ^P2 


1.45 


> 0.06 


Baldo et al., 1998, CD-Bonn potential 


lambda ^So 


0.24 


0.44-0.63 


Balberg & Barnea, 1998 



according to the quark model, the binding energy of hy- 
perons, ABaa = .Baa — 2i?A « 4 — 5 MeV, is somewhat 
less than the corresponding binding energy of nucleons, 
ABnn ~ 6 — 7 MeV. Due to these similarities, A hyperon 
pairing is expected in high density matter in analogy with 
the case of nucleons. Here, we use the A hyperon ^Sq gaps 
estimated by Balberg & Barnea (1998). The impact of A 
pairing on the thermal evolution of neutron stars has re- 
cently been investigated in Schaab et al. (1998a). 

After the onset of the superfluid state of the re- 
spective constituents of stellar matter the neutrino emis- 
sivities, the thermal conductivity, and the heat ca- 
pacity are suppressed by factors which behave like 
exp{-aTc/T) for T <C T^, where denotes the 
critical temperature and a is some constant of or- 
der unity (MaxweU, 1979, Levenfish & Yakovlev, 1994, 
Gnedin & Yakovlev, 1995, Yakovlev & Levenfish, 1995). 
Wc refer to Levenfish & Yakovlev (1994) for fitting for- 
mulas of TZsi, that are valid over the whole tempera- 
ture range T < for both types of pairing state. The 



TO = 0. In the case of m = 2 this factor behaves poly- 
nomially (Anderson & Morel, 1961, Muzikar et al., 1980, 
Levenfish & Yakovlev, 1994). The effects of to = 2 pair- 
ing on the cooling behaviour is discussed in Schaab ct al. 
(1998b) where considerable modifications were found. Due 
to the uncertainty in of the ground state quantum num- 
bers of the neutron condensate wc shall adopt here the 
more conventional case of the m = nodeless pairing. 

5.4- Photosphere 

While the heat conduction in the interior of the star is 
treated dynamically in our simulations the actual surface 
temperatures are obtained by attaching an outer envelope 
with /9 < Pm at the boundary to the 'core'. The validity of 
this approximation has been discussed in Gudmundsson 
et al. (1983). To obtain the surface temperatures we use 
the envelope calculations by Potekhin et al. (1997) whose 
non-magnetic photospheres provide a smaller gradient for 
low effective surface temperatures T^s < 3 x lO^K than 
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electron conductivity in the non-degenerate region. This 
yields a slower cooling in the photon cooling era for stars 
older than > lO^'yr. 

Even a rather small amount of accreted matter of say 
AM 10~^^Mq, substantially reduces the thermal insu- 
lation of the photosphere. The effect is highest for AM ~ 
IO^^Mq; where the accreted hydrogen is partly burned 
into helium and carbon. Additionally accreted matter will 
be converted into iron, leaving the thermal insulation 
of the star practically unchanged (Potekhin ct al.. 1997). 
The effect of accreted envelopes on cooling of neutron stars 
are discussed in Sect. 7.3. 

5.5. Comparison of the efficiency of the heating 

mechanism 

In Fig. 2 the heating rates 

H = J Awr'^e^+^'^hdr (54) 

discussed in Sect. 3 and 4 are compared with the total 
photon and neutrino luminosities and L^, respectively, 
where an isothermal configuration, i.e. Te'^ = const., is 
assumed. Since the heating rates depend not only on the 
temperature but also on the time, we use the standard 
thermal evolution scenario without heating to link the 
star's temperature to its age. 

We use two values of the rotational parameter K (see 
Eq. (5)). About the first value, K = 10~^^ s (left panel in 
Fig. 2), the observed values of most of the pulsars scat- 
ter, particularly of the five pulsars in class A. However 
there are four pulsars (see Table 3) which are rather old, 
T > 6 X lO'' yr, and rotate very fast, P < 6 ms. These 
pulsars are supposed to be recycled by accretion of mat- 
ter from a companion star (see, e.g., Lorimer, 1996). Since 
their rotational parameter K ~ 10~^^ s is quite different 
from the one of all the other pulsars, we shall study them 
separately (see Sect. 7.4 and right panel of Fig. 2). 

For the value K = 10~^^ s the rate of only three 
heating mechanism, EB- and PVB-pinning and corotat- 
ing core, exceed both, the photon and the neutrino lumi- 
nosities. We expect (and obtain) therefore only for these 
three heating mechanism a considerable effect (see Sect. 
7). All these heating mechanism affect both the late and 
the early cooling. This is in contrast to the case, where 
K = 10~^^ s is employed. Since the time derivative of the 
angular velocity Q is small in the latter case, only the late 
cooling is affected. Besides the three heating mechanisms 
quoted above, the other heating mechanisms result in a 
heating rate which is larger than the photon luminosity, 
except the heating from a corotating crust. The rate of 
chemical heating of the core depends also on the devia- 
tion from chemical equilibrium and can be compared with 
the other processes only in detailed cooling simulations. 



6. Observed data 

Only a few of the known pulsars have been detected by 
the X-ray observatories Einstein, EXOSAT, ROSAT, and 
ASCA during the last two decades. We use a sample of 27 
pulsars to compare the theoretical cooling curves with the 
observed data (sec Table 3). Both the timing characteris- 
tics (i.e. O and Q) and at least an upper limit on the effec- 
tive surface temperature as measured at infinity are 
known for these pulsars. Table 3 summarises these data. 
The effective surface temperatures are specified together 
with their 2a error range. 

The ages of all pulsars except PSRs 0531+21 (Crab), 
0833-45 (Vela) and 0002-1-62 are estimated via their spin- 
down age, T = P/2P. This relation requires that both the 
moment of inertia and the magnetic surface field are con- 
stant, and that the braking index n is equal to its canoni- 
cal value of 3, which corresponds to a spin-down by emis- 
sion of pure magnetic dipole radiation (see Sect. 2.2). Of 
course, the spin-down age is a crude approximation to the 
true age of a neutron star, which can deviate from this 
value by a factor as large as ^ 3, as discussed before. 

The situation is different for the three pulsars quoted 
above, where ages are known from different sources, that 
is: the age of the Crab pulsar is known from the chron- 
icles, the age of Vela was recently determined via the 
proper motion of the pulsar with respect to the super- 
nova remnant by Aschenbach et al. (1995), and the ap- 
proximate age of PSR 0002-1-62 is given by an estimate 
of the age of the related supernova remnant G 117.7-1-06 
(Hailey & Craig, 1995). 

The information obtained from the X-ray observations 
is not always sufficient to extract the effective surface tem- 
perature of the corresponding neutron star. The sample is 
therefore divided into four categories labelled A through 
D (Ogelman, 1995, Schaab et al., 1996): 

Category D: The four pulsars have not been detected in 
the soft X-ray range so far. By considering the instru- 
mental sensitivity an upper limit for the surface tem- 
perature could be set. These pulsars are marked with 
white triangles in the figures below. 

Category C: The detections of ten pulsars contain too 
few photons for spectral fits. The surface temperatures 
were obtained by using the totally detected photon 
flux. These pulsars are marked with black triangles. 

Category B: The spectra of eight pulsars, including the 
Crab pulsar 0531-1-21, can only be fitted by a power- 
law-type spectrum or by a blackbody spectrum with 
very high effective temperatures and effective areas 
much smaller than a neutron star surface. Presum- 
ably, their X-ray emission is dominated by magneto- 
spheric emission. Therefore, the temperatures, deter- 
mined from the spectral fits, are probably higher than 
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log(T) [K] log(T) [K] 



Fig. 2. Heating rate H as functions of internal temperature for A: EB-pinning, B: PVB-pinning, C: corotating crust, D: 
corotating core, E: crust cracking, F: chemical heating of the crust. The left panel corresponds to K = 10~^® s, the right one to 
K = 10^^^ s. The photon luminosity, as well as the neutrino luminosity of a standard cooling and an enhanced cooling model 
are also shown. 



Category A: Finally, there are five pulsars, 0833-45 
(Vela), 0656+14, 0002+62, 0630+18 (Geminga), and 
1055-52, which allow two-component spectral fits. The 
softer blackbody component is believed to correspond 
to the actual surface emission of the neutron star, while 
the harder blackbody (or power-law) component may 
be due to magnetospheric emission. These pulsars are 
marked with errorbars. 

The obtained effective surface temperature of pulsars 
of category A depend crucially on whether or not a hy- 
drogen atmosphere is assumed. PSR 0833-45 (Vela) is a 
specific example (cf. Table 3). Generally a spectral fit 
with a magnetic hydrogen atmosphere yields a substan- 
tially lower effective surface temperature than a black- 
body fit. These models of hydrogen atmospheres however 
do not yet take into account the presence of neutral atoms 
(Potekhin ct al.. 1997). Their effect should yield a spec- 
trum that is more similar to the blackbody spectrum, 
as indicated by the preliminary estimates by Shibanov et 



observations in the near future (Pavlov et al., 1996). We 
shall use the effective surface temperatures obtained by 
fitting the spectra with both the blackbody and the mag- 
netic hydrogen model atmosphere for the three pulsars of 
category A, for which both fits were accomplished. 

7. Results and discussion 
7.1. Standard cooling 

Figure 3 shows the surface temperature as measured by 
a distant observer as a function of the star's age for the 
slowly cooling neutron star models with and without 
internal heating. These models are based on the non- 
relativistic equation of state UV14+UVII which account 
only for chemically equilibrated nucleons and leptons (cf. 
Sect. 5.1). We choose the canonical value M = I.4M0 for 
the gravitational mass of a neutron star. Note that the 
relatively low surface temperature is caused by the inclu- 
sion of the superfluid pair breaking and formation process 
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Table 3. Sample of observed data 



Pulsar 


P 


P 


log(r) 


\og{K) 


log(TS) 


Category 


Reference 




[ms] 


[10-'"ss-i] 


[yr] 


[s] 


[K] 






0531+21 


33.40 


420.96 


2.97t 


-13.9 




B 


1 


(Crab) 
















1509-58 


150.23 


1540.19 


3.19 


-12.6 


6.11 + 0.10 


B 


2 


0540-69 


50.37 


479.06 


3.22 


-13.6 


6.77+°°t 

— U.U4 


B 


3,4 


0002+62 


241.81 




~ 4t 


- -13 


6.201°:"^ 


A,bb 


5 


0833-45 


89.29 


124.68 


4.3 ± 0.3t 


-14.0 


6.24 + 0.03 


A,bb 


6 


(Vela) 










5.88 + 0.09 


A,mH 


7 


1706-44 


102.45 


93.04 


4.24 


-14.0 


6.031^8 


B 


8 


1823-13 


101.45 


74.95 


4.33 


-14.1 


6.01 ± 0.02 


C 


9 


2334+61 


495.24 


191.91 


4.61 


-13.0 


0.92 Q 


C 


10 


1916+14 


1181 


211.8 


4.95 


-12.6 


5.93 


D 


11 


1951+32 


39.53 


5.85 


5.03 


-15.7 




B 


12 


0656+14 


384.87 


55.03 


5.05 


-13.7 


5.98 ± 0.05 


A,bb 


13 












t- 79+0.06 

^0 05 


A,mH 


14 


0740-28 


167 


16.8 


5.20 


-14.6 


5.93 


D 


11 


1822-09 


769 


52.39 


5.37 


-13.4 


5.78 


D 


11 


0114+58 


101 


5.84 


5.44 


-15.2 


5.98 ± 0.03 


C 


11 


1259-63 


47.76 


2.27 


5.52 


-16.0 


5.88 


C 


15 


0630+18 


237.09 


10.97 


5.53 


-14.6 


c, 7fj+0.04 


A,bb 


16 


(Geminga) 










tr 42+0.12 


A,mH 


17 


1055-52 


197.10 


5.83 


5.73 


-14.9 


5.90tn^ 


A,bb 


18 


0355+54 


156.38 


4.39 


5.75 


-15.2 


5.98 ± 0.04 


C 


19 


0538+28 


143.15 


3.66 


5.79 


-15.3 


5.83 


C 


20 


1929+10 


226.51 


1.16 


6.49 


-15.6 


5.52 


B 


21 


1642-03 


388 


1.77 


6.54 


-15.2 


6.01 + 0.03 


C 


11 


0950+08 


253.06 


0.23 


7.24 


-16.3 




B 


22 


0031-07 


943 


0.40 


7.56 


-15.4 


5.57 


D 


11 


0751+18 


3.47 


7.9 X 10"^ 


7.83 


-20.6 


5.66 


C 


23 


0218+42 


2.32 


8.0 X 10"^ 


8.66 


-21.7 


5.78 


c 


24 


1957+20 


1.60 


1.7 X 10"^ 


9.18 


-22.6 


5.53 


c 


25,26 


0437-47 


5.75 


3.8 X lO"'' 


9.20 


-21.7 


5.94 ± 0.03 


B 


27 



The entries arc; rotation period P, spin-down age t = P/2P, K = PP (see Eq. (5)), effective surface tem- 
perature as measured at infinity T^. The four categories A to D are explained in the text, bb and mH re- 
fer to blackbody and magnetic hydrogen atmosphere fits, respectively. ]\ estimated true age instead of spin- 
down age (see text). References: 1: Becker & Aschenbach, 1995, 2: Seward et al., 1983k, 3: Finley et al., 1993, 
4: Boydetal., 1995, 5: Hailey & Craig, 1995, Fig. 2 with lower limit on A^'h, 6: Ogelman, 1995, Table III, 
7: Page et al., 1996, Fig. 1, 8: Becker et al., 1995, 9: Finley & Ogelman, 1993, 10: Becker et al., 1996a, 11: 
Slane & Lloyd, 1995, 12: Safi-Harb & Ogelman, 1995, 13: Possenti et al., 1996, 14: Anderson et al., 1993, 15: 
Becker et al., 1996b, 16: Halpern & Wang, 1997, PSPC+SISO, 17: Meyer et al., 1994, Fig. 2a with Br2 = 
1.18, 18: Greiveldinger et al., 1996, Table 2, 19: Slane, 1994, 20: Sun et al., 1996, 21: Yancopoulos et al., 1994, 
22: Manning & Willmore, 1994, 23: Becker et al., 1996c, 24: Verbunt et al., 1996, 25: Kulkarni et al., 1992, 26: 
Fruchter et al., 1992, 27: Zavlin & Pavlov, 1998. 



in the earlier investigations (e.g. Umeda et al., 1994; 
Page, 1995a; Schaab et al., 1996). The standard coohng 
model without internal heating is consistent with most of 
the observations. A discrepancy occurs however for the 
effective surface temperature of the Vela pulsar 0833-45 
derived for a blackbody fit, which tends to lie considerably 
above all cooling tracks. The effective radius obtained for 
the blackbody fit, i.e. i?eff ~ 6 km, is considerably smaller 
than the canonical value of 10km, generally obtained 



reliable and in agreement with the standard cooling 
scenario without internal heating. 

The effect of heating of the crust due to thermal creep 
of pinned vortices is shown by the dotted and by the 
short dashed curves in Fig. 3 for the EB- and PVB- 
pinning model, respectively (see Sect. 3.1). The heating 
rate depends on (sec Eq. [15]) and thus on the in- 
verse of the magnetic dipole moment. The employed value 
K = lO-^S (see Eq. [5]) is in the range lO-^^^-lO"^^ s 
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Fig. 3. Standard cooling with and without internal heating. Fig. 4. Enhanced cooling of neutron star models constructed 
Parameters axe: UV14+UVII equation of state, M = IAMq, for the RHF8 equation of state. Other parameters are as in 



and K = 10 s. See Table 3 for the observed data. 



Fig. 3. 



dressed in Sect. 7.4. The effective surface tcmpciraturc is 
increased in both models in a similar way, though the ef- 
fect of the EB-pinning model is larger because of the larger 
pinning energy in the inner crust. The heating of the crust 
can easily be recognised by the increase of the thermal 
diffusion time, i.e. by the amount of time needed for the 
cooling wave to reach the surface (Lattimer et al., 1994). 
After this time the surface temperature drops significantly. 
It is increased by a factor of 50 (~ 10) in the case of 
EB- (PVB-)pinning. This characteristic feature will be- 
come very interesting if the thermal spectrum of a young 
pulsar T < 100 yr should be detected. 

For T > 10^ yr the slope of the cooling track is in- 
creased by the heating of the crust. The observed up- 
per limit of the effective surface temperature of PSR 
0950-1-08 seems to rule out the EB-pinning model. Al- 
though consideration of a magnetic photosphere would 
decrease the slope again, compensation of the heating ef- 
fect would require an unphysically strong magnetic field 
(Van Riper, 1988). For middle aged pulsars, 10^ yr < r < 
10^ yr, heating due the PVB-pinning model has only a 
small effect compared to the other uncertainties of the 
input parameters. 

The other three heating processes possible in the crust, 



are too small in comparison with the neutrino and photon 
emission rates (for the value of K used here). In the case 
of crust cracking our results are in contrast to the ones 
reported by Cheng et al. (1992), where a significant effect 
especially during the late cooling epoch has been obtained. 
This deviation can be traced back to different expressions 
used for the heating rates (see Sect. 4.1) 

Electron-vortex scattering in the interior of the star 
affects the stellar cooling at times in the range 10^ < t < 
10^ yr (see Fig. 3). During this period the effect is com- 
parable to the effect of EB-pinning in the crust. Heating 
due to vortex decay at the crust-core boundary yields an 
increase of the surface temperature only for rather old pul- 
sars (r > 10'' yr). It can not be constrained by any present 
pulsar observation, however. 

7.2. Enhanced cooling 

Figure 4 shows the cooling behaviour of enhanced cooling 
models which are based on the RIIF8 equation of state. 

The surface temperature drops by a factor of ^ 6 when 
the cooling wave has reached the surface. The fast direct 
Urea processes are only suppressed below the critical tem- 
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log(x/yr) 

Fig. 5. Standard cooling with fully accreted envelope. Other 
parameters are as in Fig. 3. 



The effect of heating on the cooling behaviour is similar 
as in the case of standard cooling. The EB-pinning model 
shows again the largest effect on the surface temperature. 

With the only exception of PSR 0630+18 (Geminga), 
the effective surface temperatures of the pulsars of cat- 
egory A are clearly not consistent with the enhanced 
cooling models, no matter whether one includes heat- 
ing processes or not. This discrepancy can be re- 
duced either by enhancing the gap energies (Page, 1995a, 
Schaab et al., 1996) or by assuming an accreted envelop. 

7.5. Accretion 

In Figs. 5 and 6, we study the effect of a fully accreted 
envelope, which means that the accreted mass exceeds 
- IO-'^Mq (sec Potckhin et al., 1997 and Sect. 5.4). Due 
to the envelope consisting of light atoms such as hydro- 
gen, helium and carbon, the temperature gradient of the 
photosphere is substantially reduced. This yields a high 
effective surface temperature in the neutrino cooling era, 
during which the energy loss rate depends on the internal 
temperature; whereas it yields small surface temperatures 
in the photon cooling era where the loss rate depends on 
the surface temperature itself. 




log('c/yr) 

Fig. 6. Enhanced cooling with fully accreted envelope. Other 
parameters are as in Fig. 4. 



class A (see Fig. 6). This is especially true for the surface 
temperatures obtained within the magnetic hydrogen at- 
mosphere model. The enhanced cooling model with heat- 
ing due PVB-pinning fits almost perfectly these surface 
temperatures. 

7.4- Millisecond pulsars 

So far we have investigated cooling scenarios associated 
with a large rotational parameter, i.e. K = 10~^^ s (see 
Eq. [5]). This value is consistent with the bulk fraction 
of observed pulsars (cf. Table 3). The four oldest pulsars, 
PSR 0751-hl8, 0218-^42, 1957-^20, and 0437-47, however 
have much smaller values around K ~ 10~^^ s. As it was 
shown in Sect. 5.5, this causes considerably different heat- 
ing rates. In particular, the mechanism of chemical heating 
and crust cracking then becomes efficient enough to no- 
ticeably influence the thermal evolution. In Fig. 7 we show 
the results for the various heating mechanism in compar- 
ison with models without heating. Since the late cooling 
behaviour (for r > 10'' yr) does not depend on the cooling 
scenario, we consider here only the standard cooling sce- 
nario. The model with chemical heating of the core shows a 
sudden rise of the surface temperature at r 10^ yr. This 
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log(T/yr) 



Fig. 7. Same as Fig. 3 but with various heating mechanism and 
K = 10"^^ s. The curves correspond to the models without 
internal heating (labeld with A), with EB-pinning (B), PVB- 
pinning (C), crust cracking (D), corotating core (E), chemical 
heating of the crust (F), and chemical heating of the core (G), 
respectively. 

idity if the deviation from chemical cquihbrium is smaller 
than the gap energy (s. Sect. 4.2 and Reisenegger, 1997). 

All cooling tracks in Fig. 7 are consistent with the up- 
per temperature limits of the three oldest millisecond pul- 
sars. Only the observation of PSR 0751+18 seems to ex- 
clude the EB-pinning model. The situation changes how- 
ever if one assumes that the actual temperatures are not 
too far from the upper limits. Then only the PVB- and 
EB-pinning models, as well as the chemical heating of the 
superfluid core yield sufficiently high surface temperature. 
A future temperature determination, which will set a firm 
lower limit on the temperature for one of these pulsars 
could therefore rule out at least the corotation and the 
crust cracking models as the only heat sources in millisec- 
ond pulsars. 

7.5. Two dimensional simulations 

In this section wc relax the simplification that the heat- 
ing rates are angle independent and extend our one- 
dimensional calculations of the previous sections to the 
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Fig. 8. Comparison of the one and two dimensional cooling 
simulations for the UVi4+UVII-model with inclusion of EB- 
pinning. 



tivity is very high in transverse direction. This is indeed a 
fairly good approximation, as we shall see next by compar- 
ing the results of one dimensional calculation with those 
of fully two-dimensional ones. The two-dimensional code 
was developed recently by Schaab & Weigel (1998) (s.a. 
Schaab, 1999). 

Fig. 8 shows the obtained cooling tracks for standard 
cooling computed for the UVi4+UVII-model with inclu- 
sion of the EB-pinning model. Heating leads to a large 
temperature difference (up to 50 %) between the pole 
{9 = 0) and the equator (9 = tt/2) for a pulsar whose 
age lies in the range 50 < t < 10^ yr. Moreover the tem- 
perature drop at the equator is delayed by the heating 
with respect to the model without heating and with re- 
spect to the drop of the polar temperature. This is the 
result of the /i oc sin dependence for the vortex creep 
model. 

The surface temperature obtained with the one dimen- 
sional code can be compared with the effective tempera- 
ture which is defined by (Page, 1995b) 

. .1 ^ 1/4 
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where s = sin S. The emission angle 6 is related to the 
colatitude 9 by the equation 



(56) 

where M and R are the neutron star's mass and radius, 
respectively. In flat spacetime, equation (56) simplifies to 

6 = e. 

The effective temperature Tcs and the surface temper- 
ature obtained with the one dimensional code are almost 
identical (see Fig. 8). This means that the luminosity ob- 
tained with both codes are also identical. The results of 
the one dimensional codes are thus acceptable for a com- 
parison with the observed data of categories B-D, where 
the upper bounds on the surface temperature are derived 
from the respective luminosity. If one uses spectral infor- 
mation as for the pulsars in category A, the observations 
should be compared with the spectra derived from the 
surface temperature distribution Ts{9) in the framework 
of two dimensional simulations. Nevertheless, we expect 
that the conclusions remain unchanged if we use the effec- 
tive temperature also for a comparison with the observed 
data of category A. 

8. Conclusions 

In this work, we have carried out a comparative analysis 
of the impact of different competing heating processes on 
the cooling behaviour of neutron stars. In general we find 
that the internal heating yields significantly enhanced sur- 
face temperatures of middle aged and old pulsars, as one 
would expect. However, the effectiveness of these heating 
processes varies significantly from one process to another, 
and depends sensitively on the value of the rotational pa- 
rameter K. 

We studied models with two different rotational pa- 
rameters. The first value was chosen to amount K ~ 
10~^^ s, which is supported by pulsars observations. For 
this K value, the heating due to thermal creep of pinned 
vortices and motion of proton vortices in the interior of 
the star yield considerable enhancements of the surface 
temperatures of middle aged pulsars with respect to the 
models which ignore heating. This leads to closer agree- 
ment with the observed data in the case of enhanced cool- 
ing, which is even more improved for the case of the fully 
accreted envelopes (s.a. Page, 1996). 

The observed upper temperature limit for PSR 
09504-08 seems to rule out the strong pinning models, 
whereas the weaker pinning models are consistent with 
the observations (the pinning models differ by the value of 
the pinning energy in the inner crust). The other heating 
processes - from chemical heating to dissipative motion 



due to crust cracking is in contrast to the one obtained by 
Cheng et al. (1992), who find a large effect especially on 
the late cooling stages. The difference arises from different 
expressions adopted for the heating rates. 

The four oldest pulsars of our sample rotate very fast 
with periods of a few milliseconds. Hence the rotational 
parameter for millisecond pulsars is smaller than that for 
the bulk fraction of the population by 7 orders of mag- 
nitude {K = 10~^^ s). For this value only the late cool- 
ing stages arc; aff(!cted by heating. We find that the heat- 
ing due to thermal creep of pinned vortices and chemical 
heating of the core has the largest impact on the surface 
temperatures of millisecond pulsars. Again, the strong pin- 
ning models leads to surface temperatures which are larger 
than one of the upper limits deduced for the millisecond 
pulsars. 

As far as we know, the published cooling simulations 
of neutron stars which account for internal heating were 
performed in one spatial dimension. Such a treatment is 
only possible if the angle dependent heating rates are av- 
eraged over spherical shells. By comparing the outcome 
of one-dimensional simulations with fully two-dimensional 
simulations, based on a recently developed numerical code 
(Schaab & Weigel, 1998), we could confirm the validity of 
this procedure. 

The knowledge of the input parameters for cooling sim- 
ulations such as the composition, the superfluid gap en- 
ergies, and the possible existence of an accreted envelope, 
etc. is far form being complete. The number of param- 
eters involved in the cooling simulations increases when 
the internal heating processes are taken into accoimt. In 
principal, one could expect to derive some of the param- 
eters associated with these phenomena by means of com- 
paring the theoretical cooling models with the body of 
observed data. This attempt needs to be supplemented 
with the analysis of terrestical experiments, other astro- 
physical observations, as well as future theoretical stud- 
ies. Nevertheless, the observed data already constrain the 
microphysical input, as, for instance, for the scenaria of 
heating via vortex creep. 

Tables containing detailed references to the in- 
gredients used in the simulations, the observational 
data, and the resulting cooling tracks can be found 
on the Web: http://www.physik.uni-muenchen.de/sektion 
/suessmann/astro/cool. 
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